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1. Introduction. In many gasdynamical problems one tries to calculate the steady state 
solution by solving the corresponding time dependant problem. One hopes that for t — ► oo the 
solution converges to a unique steady state. Recently, M. D. Salas, S. Abarbanel and D. Gottlieb 
[l]considered the initial-boundary value problem 


u t + ^(u 2 ) x = f(x), t> 0, 0<£<7T, 

u(£,0) = g(x ). 


(LI) 


They used 

f(x) = sinxcosrr, g(x)=bsmx, 0 < 6, 

and showed that the solution u(x, £) of the above problem converges to a steady state v(:r), as 
t — ► oo, but that depends on the initial data. 

In this paper we consider the viscous problem 

u t + ^(u 2 )z = ™ xx +f(x) } t> 0, 0<£<1, 6>0, (1.2a) 

with initial and boundary conditions 


u{x,0) =g{x), 
u(0, i) = a, ii(l, t) = 6, 


and the corresponding steady state problem 

0 < x < 1, e > 0, 
y(l) =6. 


2 {y 2 )x = eyxx + f(x), 


3/(0) = a, 


(1.26) 


(1.3) 


For simplicity we restrict ourselves to two cases: 


1) a > 0 > 6, a > —6, f(x) = 0, 

2) a = b = 0, f is such that there exists an a with 0 < a < 1 such that f(x) > 0 for 0 < x < a, 

f(x) < 0 for a < x < 1, /( 0) = /( 1) = 0, f x ( 0) > /o > 0 and /*(!) > /q. 


We will show that (1.3) has a unique solution and discuss the properties of y(x). We shall 
also show that in all cases we consider, the limit of y(x) as e — ► 0 exists. Thus, if 

lim ti(x, t) = y(x) 

t— ►oo 
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exists, we obtain a unique steady state solution of the inviscid equation (1.1) if we first let 
t — *■ oo and then e — ► 0. This is in contrast to the procedure in ( 1], where the two limit 
procedures are taken in the reverse order. 

We shall prove that the eigenvalues of the eigenvalue problem 

M> = ~{w>)x +e<Pxx, £>(0) = <p( 1) = 0, (1.4) 

are all negative. Therefore, the solution of (1.2) converges to the solution of (1.3) provided 
u(x,0) = g(x) is sufficiently close to y(x). In another paper we shall prove that u(x, t) converges 
to y(x) as t — ► oo for arbitrary initial data. The speed of convergence is determined by the 
eigenvalues, Ay, of (1.4). We shall show that the eigenvalue distribution depends on /(x) and 
on a, b in the following way: 

There is a constant c > 0 which does not depend on e such that 


(1) if a > — 6, / = 0 then 0 > -c/e > Ai > A 2 > • • • 

(2) if a — —6, / = 0 then - Ai = 0(e~ l / c ) > 0, -c/e > A 2 > A 3 > - • • 


(3) 

if a = b = 0, 

f f{x)dx ^ 0, 
Jo 

then — c > Ai > \ 2 > • • • 

(4) 

if a = b = 0, 

[ f(x)dx = 0, 
Jo 

then — Ai = 0(e _1 / ff ) > 0, 


(1.5) 


We expect a reasonable speed of convergence in the first and third case, while in the second 
and fourth case the speed should be extremely slow due to the eigenvalue — Ai = 0(e -1 / ff ). This 
is confirmed by numerical experiments. We see that at first u(x, t) quite rapidly approches the 
same limit as the inviscid equation (1.1), which consists of solutions of the stationary equation 


^(« 2 )* = /(*) 


connected by a shock. Once the viscous shock has been formed, the solution of (1.2) becomes 
quasi-stationary and the shock creeps extremely slowly to the “right” position. We can ex- 
plain the behavior, because by linearizing around the quasistationary solution we find that the 
eigenvalues of the corresponding eigenvalue problem have a similar distribution as earlier. 

If -Ai = 0(e~ l / e ) then the speed of convergence is so slow that the above method to 
calculate the steady state is impractical, see figures (1) and (3). However, we can use the same 
technique as Hafez, Parlette and Salas in [2] to speed up the convergence. See figures (2) and 

(4). 

Unfortunatly, not only the speed of convergence but also the condition number of the 
stationary problem deteriorates. We have to calculate with 0(e 1 / e ) decimals to obtain correct 


2 



results. To avoid an excessive number of decimals we have used a quite large e in our numerical 
calculations. 

The situation becomes much better in a two dimensional case, which we discuss in the last 
section. Now there is a whole sequence of eigenvalues 

~H!j =0(j 2 e), j = 1,2,..., 

close to zero. However, they are only algebraically and not exponentially close to zero. We 
indicate how to modify the procedure to accellerate the speed of convergence. 

We believe that the viscous model (1.2) better explains what happens in actual calculations 
than the inviscid equation (1.1). Practically all numerical methods have some viscosity built in. 
Also, from a physical point of view, the solution we are interested in is the limit of solutions of 
a viscous equation. 

Finally we want to point out that the appearance of small eigenvalues has also been 
observed by D. Brown, W. Kath, H. O. Kreiss and W. Henshaw, M. Naughton (private com- 
munication). 


2. Uniqueness, existence and properties of the steady state solution. We start 
with uniqueness, which can be proven by standard techniques. 

Lemma 2.1. If the steady equation (1.3) has a solution, then it is unique. 

Proof. Let u, v be two solutions. Then w = u — v is the solution of 

^(pw) x = ew XXl p = u + v, w(0) = u;(l) = 0. (2.1) 

If w ^ 0 then the zeros of w are isolated. Let x with 0 < x < 1 be the first zero to the right 
of x = 0. Without restriction we can assume that w > 0 for 0 < x < x, i.e. ^(0) > 0 and 
w x {x) < 0. Integration of (2.1) gives us 


-e(K(*)l + K(°)l) = e Klo = ^Wo = 0. 


Thus u^O) = w x (jc) = 0. We can consider (2.1) as an initial value problem with initial data 
iu(0) = u>z(0) = 0 whose solution is w(x) = 0, and the lemma is proved. 

We shall now discuss the properties of the solution. Let us start with the case f(x) = 0, 
a > 0 > b, a > —b. Integrating (1.3) gives us 


eyx = ^y 2 -c, 0 < a; < 1, 

y(o) = a. 


( 2 . 2 ) 
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The constant c has to be determined so that y(l) = b. We necessarily have c = d 2 / 2 > a 2 / 2, 
because with c < a 2 / 2, y x > 0 for all x, and y(l) = 6 cannot be satisfied. We can solve 
equation (2.2) explicitly. This is done by writing (2.2) in the form 



a 0 


i.e. 


( 


a + d 

a — d 


)( 


y{x)+d’ 


Therefore y(l) = 6 implies d = a + 0(e -1 / ff ), and 


y(x) - a 


1 _ Te -a(l-x)fe 

1 + 


with 


a — 6 
a + 6 


(2.3) 


Away from the boundary layer at x=l we have y(x) = a + 0(e~ a ^~ x ^ c ). Thus, for e — * 0, y(a;) 
converges to a for 0 < x < 1. 

If a = - b we consider (2.2) on the interval 0 < x < i, with boundary conditions y(0) = 
a, y(^) = 0 and obtain a solution yi(:r) of the form (2.3). The solution on the whole interval 
is given by 


y(x) = 


\ yi(s), 

\ — 2/i(l 


x), 


if 0 < x < 1 
if 5 < * < 1. 


In figures (9) and (10) we have plotted y(x) for two different sets of boundary values. 

Consider case 2, where / only vanishes at x = 0, a, 1 and a = b = 0. Without restrictions 
we can assume that 

l 

J f(x) > 0. (2.4) 

0 


If this is not true, we transform the problem by introducing new variables, 

x = 1 - x, / = -/, y = -y. 


The new problem satisfies (2.4). 

Lemma 2.2. Let y(x) be the solution of (1.3), F(x) = f(£)d£ and 

Then 

y*{ 1) < Vx{0) < K u Ki = ma^c {\h x (x)\} + |h,(0)|. 


h(x) = y/2F(x). 
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Proof. Integration of (1.3) gives 


e{y x - y*(0)) = |y 2 - F, 


(2.5) 


y(o) = o, 


where y x (0) is determined by y(l) = 0. If u = y - h, then u is the solution of 

tlx = yj(0) hx ”i” £ uh d - 2 £ ^ > 


u( 0 ) = 0 . 


Assume that y 2 (0) > K\. It follows that y x {0)-h x (x) is positive and thus u and u x are positive 
for all x > 0. In particular u(l) > 0 and y(l) = u(l) + h( 1) > 0, which contradicts y(l) = 0. 
Thus yz(0) < K\. Also 

ey a; (l) = ey x (0) - F( 1) < ey x { 0). 

This proves the lemma. 

Lemma 2.3. Let y(z) be the solution of (1.3) and let e be sufficiently small. If F(l) > 0 
then y(x) > 0 for 0 < x < 1 and y(x) has exactly one maximum. If F(l) = 0 then there exists 
an x with 0 < x < 1 such that y(i) > 0 for 0 < x < x, and y(x) < 0 for x < x < 1. Also y(x) 
has exactly one minimum and one maximum. In both cases |y(x)| < max|F(z)|. 

Proof. At extrema y x = 0 and 


Vxx = -e l f = 


<0 for 0 < x < a, 
= 0 for x = a, 

> 0 for a < x < 1. 


( 2 . 6 ) 


Thus y cannot have a minimum to the left of a maximum. Since y(0) = y(l) = 0 there are only 
three possibilities, namely 


y > 0 

for 

0 < X < 1, 

y has exactly one maximum, 

(2.7a) 

y < 0 

for 

0 < x < 1, 

y has exactly one minimum, 

(2.76) 

y > 0 
y < 0 

for 

for 

0 < x < x, 

X < X < 1, 

0 < x < 1, 

y has exactly one maximum and one minimum. 

(2.7c) 
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We shall prove that if F(l) >0 then (b) and (c) are not possible, and that if F(l) = 0 
then (a) and (b) are not possible. 

Let F( 1) >0. Suppose (2.7b) holds. Then 


y*(o) < 0, y x (l)>0. 


By (2.5) 


0 < e(j/i(l) — !/*(0)) = —F(l) <0. 


( 2 . 8 ) 


This is a contradiction, so (2.7b) cannot hold. Now supposse (2.7c) is valid. Then y-^O) > 0 
and by (2.8) 

y*(0) > e -1 F(l). 


If e is small enough this is impossible by lemma 2.2. 

Let F( 1) = 0. Assume that (2.7a) or (2.7b) are valid. By (2.8) y*(0) = y*(l), which is 
only possible if y*(0) = y^l) = 0. Differentiating (1.3) gives us 


£ !/in — llllxz + (yx)~ fx- 


(2.9) 


Thus 


y(o) = y*(o) = y xx {o) = o, y IM (o) < o, 

y(f) = J/a(f) = 2/is(l) = o, y xxx {l) < 0. 

This implies that y must change sign at least once, which contradicts the assumption, and 
therefore (2.7c) must hold. 

It remains to show that |y(x)| is bounded by max|F(x)|. Since y(0) = y(l) = 0, the 

maximum absolute value of y is found at a local extrema, where y x = 0. Thus, from (2.5) it 

follows that 

M*)! ^ “2* l-Ffc) “«y*(0)| < max |F(x)|. 


This finishes the proof. 

We can use the usual singular perturbation methods to discuss the behavior of the solution 
in detail, see for ex. [3]. 

Theorem 2.1. Let F(l) > 0, assume that (1.3) has a solution and that e is sufficiently 
small. Then y(x ) has a boundaiy layer at x = 1. For 1 - 0(e|log(e)|) < x < 1, y(x) is close to 
w(x) which is the solution of 


ew x = -w 2 - F( 1), 


—oo < x < 1, w(l)=0. 


In any interval 0<x o <x<l — 0(e\ log(e)|) 


y(x) = /i(x) + sui (x, e), h(x) = F(x) =: xg(x), 


(2.9) 


( 2 . 10 ) 
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where u\ and its derivatives are bounded independantly of e. For 0 < x < xo < a we have 

y(x) = h(x) + eu(x), x = x/y/e ) (2.11) 

where u and the derivatives d v u/dx v are bounded independantly of e. Thus, for e — * 0, y(x) 
converges to h(x) for 0 < x < 1. 

Proof. We indicate only the proof of (2.11). In the proof we shall use l\,h and 1 to 
denote the intervals 0 < x < 1, 1 < x < Xo/y/e and 0 < x < xo/y/e, respectively. We shall also 
use 

11/11/ :=max|/(x)|, 


where I is an interval. 

We introduce a new variable in (1.3), 

y(x) = h(x) + eu(x/y/e). 


This gives us 

v>xx — {%g{ x ) + y/eu) u z — h x u — — h xx , 0 < x < xo/ y/e, t/(0) = 0, u(xq/ y/e) — u 0) (2.12) 

where uo = ui(a:o,e) is bounded independantly of e. From xq < a and the assumtion f x ( 0) > 
/ 0 > 0 it follows that h x {x) > h 0 > 0 for 0 < x < x 0 . Therefore we can use the maximum 
principle. The maximum of u is found either on the boundary or at a local extrema, where 
Ux = 0. At local extrema 

M < 1^1 < 


Thus 

||ti||/ < max(tto,o:). (2.13) 

Next we want to estimate ||tt 5 ||/. First we consider the interval I\ = [0, 1]. By (2.12) and 
(2.13) there are constants C\ and Ci such that 

\WxxWh < CllWxWh + 

It is well known, see Landau [4], that one can estimate \\uz\\r l in terms of ||u|| and ||if££||/ t , 
i.e. for every constant 8 there is a constant C(8) such that 

Thus for 8 = ±(Ci) _1 we obtain a bound for ||tia«II/n which gives us a bound for Hu^H/j. 
Especially, |ti£(l)| is bounded. 

In the remaining interval I 2 = [1,®o/v^]j we have 

F>F(V^)=ef x (0)(l+O(y/i)). 


Thus 

xg + y/eu = \PlF j y/e + y/e > \J f x { 0 ) + 0{y/e), 


7 



i.e. for sufficiently small y/e 


xg + s/eu> • 

At local extrema of u x , uj* = 0 and we have, by (2.12), 


^ “ \/]x(0)^ hXX ~ hxU \ ~ s/7Jd)^ hxX ^ h + =: P- 


Thus 

INII/* < max(|«i(l)|, |ui(-^)|, 0), 

V e 

and u~ x is bounded independantly of e in the whole interval. By differentiating (2.12) bounds 
for higher derivatives of u can be obtained. 

It is also clear that as e — ► 0, y(x) converges to h(x). This finishes the proof. 

If F(l) = 0 then the solution switches at x from \/2F + O(e) to -\Z^F + 0(e). In each 
subinterval 0 ^ x < x and x < x < l the local behavior of the solution is of the same type as 
in the first case. As s -*• 0, y(x) converges to h(x) for 0 < x < x and to ~h(x) for x < x < 1. 
In general, the position of x can only be obtained by detailed calculation. However, if f(x) is 
antisymmetric around x = | then x = \. This is the only case we consider. 

We shall now discuss the existence of a solution. For this we need two lemmata. 

Lemma 2.4. For sufficiently large e the steady state equation (1.2) has a solution, 
Proof. By integrating (1.3) twice, we can write the equation in the form 


1 x X 

y (*) = j yf y 2 (Z) d Z -nj + vxc 0 , v = i/e, 


1 «+co = o, 
0 0 

or after the change of variable y = r\y 


1 X X 

y(a) = l*l 2 f y 2 (0dt- 1 F(Z)d£ + xc 0 , 


\v 2 1 J F(Z)dZ + co = 0. 

o 0 
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For t) — 0 the above equations have a unique solution. Therefore the same is true for all 
sufficiently small p. This proves the lemma. 

Lemma 2.5. Let p(x) be a smooth function. Consider the eigenvalue problem 

\<p = ~(p<p)x + e<Pxz, <p(0) = <p(l) = 0. (2.14) 

The eigenvalues are real and negative. 

Proof. We introduce a new variable xp(x) by 


5 e-, /p(04£ 

<p{x) = e 1 xp(x), 

and obtain 

At f> = ex p xx ~ exp =: Lxp, c(x) = \p x {x) + -^(p(x)) 2 , 

L 4S (2.15) 

xp(0) = xp{\) = 0. 


(2.15) is selfadjoint and therefore the eigenvalues are real. Let <p ^ 0, A be a solution of (2.14), 
and let x be the first zero of <p to the right of x = 0. We can assume that <p > 0 for 0 < x < x. 
Thus (p x { 0) > 0 and <p x {x) < 0, and integration of (2.14) gives us 


A j <p(x)dx = e[<Px]o < 0. 
o 

It follows that A < 0. If A =0, the only possible solution of (2.14) would be (p(x ) = 0. Thus 
A < 0, which proves the lemma. 

Now we can prove 

Theorem 2.2. The equation (1.3) has a unique solution for all e > 0. 

Proof. We have already shown that (1.3) has a solution for sufficiently large e . We will 
now employ continuation in e to prove existance for all e > 0. Assume we have shown existance 
for e > 6 . We want to show that there is a solution for e — e . By lemma 2.3 the solutions of 
(1.3) are uniformly bounded for e < e < e + 1. Therefore the same is true for the first three 
derivatives. Thus we can select a sequence of solutions 

y{x,e v ), v — 1, 2, ... , = e, 


such that 


7 y(x, et/ ) = £jy(x,e), j = 0,1,2 , 
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and y(x,e) is the desired solution. Linearizing the equation around y(x,e ) gives us 

(y(x,e)6y) x = e(6y) xx + (e -e)y(x,e), 6y(0) = 6y(l) = 0. ^ 

By the previous lemma A = 0 is not an eigenvalue of the above equation and therefore we can 
solve (1.3) for all sufficiently small e — e. This proves the theorem. 


3. Speed of convergence. In this section we want to discuss the speed of convergence 
to steady state. We assume that the initial data g(x) of (1.3) are sufficiently close to the 
solution of the steady problem, so that we only have to discuss the behavior of the solutions of 
the linearized equation 


w t + {yw) x — ew xx , 0 < x < 1, t > 0, 

m(x,0) =ff(ar), (3.1) 

u)(0, t ) = u/(l, t) = 0. 


To determine the speed of convergence we study the distribution of eigenvalues of 

* <P + (y<P)x = e<p xx , .p(O) = ^(1) = 0. (3.2) 

Theorem 3.1 . The eigenvalues of (3.2) are real and negative and their distribution is 
given by (1.5). 

Proof . Lemma 2.5 tells us that the eigenvalues are real and negative. First we consider 
the case f = 0, a > -b. We write (3.2) in the selfadjoint form (2.15) with p = y. Let A = Ax 
be the largest eigenvalue. The corresponding eigenfunction V’l does not change sign, and we 
can assume that ipi > 0 for 0 < x < 1 and that max |-0i (^c)| = 1. We assume that Ai > —a 2 / 8e. 
Then there is a constant K such that c(z) + A x > 0 for 0 > x > 1 - Ke. Thus rpi is monotone 
in the interval 0 < x < 1 — Ke, and therefore tpi must have its maximum in the remaining 
interval, 1-Ke <x <1. By assumption maxt/>x(x) = 1 and therefore there must be a constant 
6 > 0 such that 1 ) < —8/e. Now consider the corresponding eigenfunction 


fy(m 

<Pi{x)=e « V’i(ar), = rp lx {l), 0 < ^(x) < ^(x). 

Integrating (3.2) gives us 


-6 > e(<p lx (l) — (p lx (0)) = \i [ cpi dx> 

Jo 


>A, I' e*'~'K yd tdx = \ l ed. 

Jo 
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Thus 


\ . ,a 2 <5 . 

A, 


and the theorem is proven for this case. 

When / ^ 0, a = 6 = 0, and / 0 l f(x)dx > 0 the corresponding estimate follows in the same 
way, since by theorem 2.1 there are constants Cb ^ 0 and K such that 

c (x) = i(/i 2 (x)+0(^) + ie _I (/i(x)+0( v /£)) 2 ) >C 0 >0 for 0 < x < 1 - Ke. 


We now consider the antisymmetric case when a — —b, f = 0 or a — 6 — 0 and f(x) is 
antisymmetric around x = We want to show that 


— Ai = 0(e l e 

We shall use the fact that for our selfadjoint eigenvalue problem (2.15) the eigenvalue with the 
smallest absolute value, Ai , satisfies 


M< 


¥4h 

Uh ’ 


for any smooth function <f>£ 0 satisfying the boundary conditions. We chose 


%*~ l / y (€)<*€ / y(0<t£ 

<f>(x) = e 1/2 — e 0 

as trial function. y(x) is antisymmetric around x = and <f>(0) = </>(!) = 0. Also 


o 1 1 

r / ,y , y *\. 2 
^ ( fe + T )e 


1/2 

f y(0<*€ 


0 


Both <{> 2 and {\e~ l y 2 + |yi) 2 are symmetric around x = |. Therefore 


1/2 


IMP = 2 /(g + f)V‘-‘ 


1/2 


ll^l 2 = 2 f e" 5 "' fo /2 ^(e^~' fo - l) 2 dx, 
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and by (2.5) and theorem 2.1 


1/2 2 

{ (g + 

0 


y2<J!MJ!_ 

1_ W " 1/2 


< C 2 e~ 2 e~ 2D l e , 


/ (e 

0 


— l) 2 cf:r 


where C > 0 , D > 0 are constants which do not depend on e . 

We shall now estimate the size of the second eigenvalue for the case with an interior 
boundary layer at x = ±. By assumption y(x) is antisymmetric around x = A. Consider the 
eigenvalue problem (3.2) on half the interval, 0 < 2 < and denote its solutions by 


cpi(x), Xi, i = l,2 

We know that p, has i - 1 sign changes, and we have already shown how the A,-’s are bounded 
away from zero. The function 


<P2i(x) = 


<pi ( 2 ) for 0 < x < | 

-<Pi(x — |) for \ < x < 1 ’ 


* = 1 , 2 ... , 


will satisfy (3.2) on the full interval, 0 ^ x ^ 1 with A — A 2 i — A,. Also changes sign 
2(t — 1) + 1 times. Thus <p- 2t is the 2 i th eigenfunction and X 2 i is the 2 i th eigenvalue. Therefore 
A 2 is bounded away from zero. This finishes the proof. 


4. Numerical results. We shall discuss difference approximations for the time depen- 
dant problem (1.2) and the eigenvalue problem (3.2). We introduce gridpoints 


( 2 , = ih, tj = jk), i = 0,1,... y = 0,l h 


_ 1 _ 

N* 


where N is a natural number and k > 0 is the time step. We also introduce gridfunctions 

«i =«(* 

We approximate (1.2) by the usual implicit method 


(7 - ekD + D-)ui +1 + hD 0 (x4 +1 ) 2 = uj'+ kf i} i = 1, 2, . . . , N - 1 


( 4 . 1 ) 
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with initial and boundary conditions 


u® — (ft, t — 1,2, ...,7V 1, 

= a, = 6, J = 1,2,. . . 

Here 

h 2 D+D-Ui = tii + i - 2 Ui + Ui_! and 2/z£>o(u,*) 2 = («*+i) 2 ~ (t«i-i) 2 

denote the usual centered difference operators. At every time step one has to solve a nonlinear 
system to determine u% +l . This is done by the iteration 

(I - ekD+D-)ri l+1) = ~kD 0 (u<P) 2 + u>+ kf it 1= 0,1... , (4.2) 


where is choosen by a predictor process. 

In all our experiments the solution of (4.1) converges to a steady state solution. However, 
the speed of convergence depends on the location of the shock. If the shock is located at the 
boundary, corresponding to the first and third case of (1.5), then the convergence to steady 
state is quite rapid. See figure (5). If on the other hand the shock is located in the interior, 
corresponding to the other cases of (1.5), the convergence is, in general, very slow. When the 
shock is formed at an early stage it is in general in the “wrong” place, depending on the initial 
data. From then on, the the shock moves slowly to the correct position. See figures (1),(3). This 
process can be considered quasi-stationary, which makes it possible to use the same convergence 
acceleration as in [2]. 


Formally we can write our iteration (4.1) as 


H{u n+1 ) = u n+1 -u n := r n . 

(4.3) 

We can linearize the realation and obtain 


(/ - L)r n+1 = r n . 

(4.4) 

In our case 


Lri = ekD+D-ri — kDo(u™ +l r t ). 

(4.5) 


This is a discretization of the right hand side of the eigenvalue problem (2.14), with p = u n . If 
the process is quasi-stationary we can consider L to be independant of n. Then we have 


r n+i = {I-L)-’r n 


and 


p-i 


u n+P = u n + J2(I-L)- j r n . 

)=0 


If the eigenvalues A,-, of L are negative the eigenvalues k,, of (/ — L) 1 satisfy |/Cj| < 1 and 

lim u n+p = u n + (/ — (/ — L)" 1 r 1 r" = u n + (/ — L~ l )r n . (4.6) 
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Instead of taking a large number of time steps we can take one large step, which we call an 
extrapolation step. We put 

u = u n + Pe, (4.7) 

where e is the solution of the equation 

Le=(L- 1) r", (4.8) 

and p is a stabilizing parameter. We choose P in such a way that H(u n + fie) has no component 
in the direction of e, i.e. 

(H(u n +/3e),e)=0, 

where {•, •) denotes the usual inner product. There are other possible choices, for example 
choose /? such that 

||i/(u" +pe) || = nun + /?e)||. 

P 

Of course (4.7) is not the steady solution we are seeking. We use the new u to restart the time 
iteration, and make a new extrapolation step once a new quasi-stationaiy state is reached. In 
our experiments we use an a priori fixed number of time steps between the extrapolation steps. 
Better strategies are under development. 

We have calculated the first eigenvalues and eigenvectors of the discrete linearized operator 
(4.5), provided u" +1 is the discrete steady state solution. The calculations show that the 
eigenvalues are negative and their distribution is of the same type as for the corresponding 
continous case. See table (1). In figures (6), (7) the first few eigenvectors are plotted. Note that 
in the case of an interior shock the first eigenvector is exponentially small away from the shock 
region. Also, we have no doubt, and it is confirmed by the calculations, that the position of the 
shock does not change the nature of the eigenvalue distribution. In fact, in the proof of theorem 
3.1, y can be replaced by any function of the same structure. 

In our case, when the shock is located in the interior, (I — L)~ l has only one eigenvalue, 
/ci , close to zero. All other eigenvalues are small. Therefore, when we have reached the quasi- 
stationary state, r n is in the direction of the eigenvector corresponding to /cj. See figure (8). 
Therefore we do not need to solve (4.8), and instead of (4.7) we use 


u = u" + pr n . (4.9) 

In figures (2), (4) we have plotted u at different time stages to show how the convergence is 
accelerated. 
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5. A twodimensional case. Consider the following problem 


u t + = e(u xx + Uyy), 0 < x < 1, 0 < t/ < 1, t>0, 

£ 

u(0, y, t) = a, u{l,y,t) = -a, a> 0, 
u(x,0, t) = u(x,l, t) = w(x), 
u{x,y,0) =g{x,y), 


(5.1) 


where W(x) is the solution of the one dimensional problem (1.3) with 6 = —a, and f(x) = 0. 
See (2.3). A steady solution of (5.1) is u(x,y) = w{x). 

The speed of convergence can be studied by analyzing the corresponding eigenvalue prob- 
lem 

H<p + (w<p) x = e(ifxx + 'Pyy)) <P = 0 on the boundary. (5.2) 

We can solve (5.2) by separation of variables. Let <p[x,y) = X(x)Y (y). Then 

(u;X) ; - eX" = \X, X(0) = X(l) = 0, (5.3a) 

Y" ---- -qY, F(0) = F(l) = 0, (5.36) 

with (i — \ — £q. We recognize (5.3a) as (3.2). Therefore — Aj = 0(e -l / e ) and —\j > 0{ 1/e), 
j = 2,3, We can solve (5.3b). The solution is 

Yj(y) = sin (j Try), qj = (jnf, j = 1,2 

There is a whole sequence of eigenvalues, /Xjy, of order O(e). The eigenfunctions corresponding 
to this sequence, (p\j, will be exponentially small away from the shock. All other eigenvalues 
will be of order 0(\/e). 

We expect that the time iteration will again lead to a quasi-stationary state, and that 
the residual will be composed of eigenfunctions corresponding to the eigenvalues of order 0{e). 
Therefore e in (4.8) will be of the same form, and we can replace all components of e away from 
the shock by zero, thus obtaining a linear system of equations of order N instead of N 2 . More 
details will be given in another paper. 
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Table 1. 

Eigenwlues of the eigem'alueproblem (3.2), y is the solution of (1.3). Three different cases were treated. 
The discretization is done according to (4.5), with N — 100 gridpoints . The eigenmlues were found using 
inverse iteration. Eigenvectors corresponding to case (1) are plotted in Ggure (6a, b). 



■■ 

u 

*3 

/(x) = sin(27rx)/2 
a = 6 = 0 
e = 0.04 

-8.64 • 10~ s 

-4.34 

-5.32 

/(x) — sin(2rrx)/2 
a = 6 = 0 
c = 0.02 

-4.62 -lO" 6 

-5.617 

-5.622 

f(x) = 0 
a= 1, 6= -1 
£ = 0.02 

-1.24- 10~ 9 

-12.8 

-13.5 
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Figure 1. Convergence in time without convergence acceleration. Numerical solutions at 
different time stages for the case e = 0.05, / = 0, a = 1, b = —1, u(x, 0) = 1 + 2(e~ 2x — 1) /(I — e“ 2 ). 
Between each curve there are 200 time steps = 40 time units. The calculation is made with time step k 
= 0.2 and N=50 grid points. 
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Figure 2. Convergence in time with convergence acceleration. Numerical solutions at different 
time stages for the same case as in figure 1. Between each curve there are 15 time steps and one 
extrapolation step. The same time step, k=0.2, and number of grid points , N=50, are used. 
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Figure 3. Convergence in time without convergence acceleration. Numerical solutions at 
different time stages for the case e = 0.04, / = - sin(7rx) cos(ttx), a = b = 0, u(x, 0) = ^sin(7rx). 
Between each curve there are 100 time steps . The calculation is made with time step k = 0.1 and N=50 
grid points. 
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Figure 4. Convergence in time with convergence acceleration. Numerical solutions at different 
time stages for the same case as in figure 3. Between each curve there are 20 time steps and one 
extrapolation step. The same time step, k=0.1, and number of grid points , N=50, are used. 
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Figure 5. Convergence when the shock is located at the boundary. Here e = 0.04, f(x) = 
- sin(7rx), u(x, 0) = ^sin(7rx), N = 50, k = 0.1. Between each curve there are 5 time steps. 
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Figure 6a. Eigenvectors. The first two eigenfunctions of problem (3.2), when y, the solution of (1.3), 
has a shock in the interior. In this case e = 0.04, f(x) = j sin(Trz) cos(7rz), a = b = 0, N = 100. 
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Figure 6b. Eigenvectors. The third and fourth eigenfunctions of problem (3.2), when y, the solution 
of (1.3), has a shock in the interior. In this cases = 0.04, f(x) = ~ sin(7rx) cos(7rx), a = b = 0, N = 
100 . 
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Figure 7. Eigenvectors. The first two eigenvectors, <pi and <p 2 , of problem (3.2), when y, the solution 
of (1.3), has a shock x = 1 . In this case e = 0.08, f(x) = - sin(7rx), a = b = 0, N = 100. 
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Figure 8. Differences between consecutive solutions at different time stages, when e = 0.04, 
/ = ^sin(7rx)cos(7rx), a = b = 0, u(x, 0) = ^sin(7rx). Between each curve there are 100 time steps . 
The calculation is made with time step k = 0.1 and N=50 grid points. 
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Figure 9. The solution of (1.2) when / = 0, a = l,b = 0ande = 0.05. 
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Figure 10. The solution of (1.2) when / = 0, a = 1 , 6 = — 1 and e = 0.05. 
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